GCC Code Coverage Report


Directory: ./
File: phys/phyaqua_mod.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 382 0.0%
Branches: 0 488 0.0%

Line Branch Exec Source
1 !
2 ! $Id: phyaqua_mod.F90 3579 2019-10-09 13:11:07Z fairhead $
3 !
4 MODULE phyaqua_mod
5 ! Routines complementaires pour la physique planetaire.
6 IMPLICIT NONE
7
8 CONTAINS
9
10 ✗ SUBROUTINE iniaqua(nlon,year_len,iflag_phys)
11
12 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13 ! Creation d'un etat initial et de conditions aux limites
14 ! (resp startphy.nc et limit.nc) pour des configurations idealisees
15 ! du modele LMDZ dans sa version terrestre.
16 ! iflag_phys est un parametre qui controle
17 ! iflag_phys = N
18 ! de 100 a 199 : aqua planetes avec SST forcees
19 ! N-100 determine le type de SSTs
20 ! de 200 a 299 : terra planetes avec Ts calcule
21
22 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23
24 USE dimphy, ONLY: klon
25 USE geometry_mod, ONLY : latitude
26 USE surface_data, ONLY: type_ocean, ok_veget
27 USE pbl_surface_mod, ONLY: pbl_surface_init
28 USE fonte_neige_mod, ONLY: fonte_neige_init
29 USE phys_state_var_mod
30 USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
31 day_ini,day_end
32 USE indice_sol_mod
33 USE nrtype, ONLY: pi
34 ! USE ioipsl
35 USE mod_phys_lmdz_para, ONLY: is_master
36 USE mod_phys_lmdz_transfert_para, ONLY: bcast
37 USE mod_grid_phy_lmdz
38 USE ioipsl_getin_p_mod, ONLY : getin_p
39 USE phys_cal_mod , ONLY: calend, year_len_phy => year_len
40 IMPLICIT NONE
41
42 include "YOMCST.h"
43 include "clesphys.h"
44 include "dimsoil.h"
45
46 INTEGER, INTENT (IN) :: nlon, year_len, iflag_phys
47 ! IM ajout latfi, lonfi
48 ! REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
49
50 INTEGER type_profil, type_aqua
51
52 ! Ajouts initialisation des surfaces
53 ✗ REAL :: run_off_lic_0(nlon)
54 ✗ REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
55 ✗ REAL :: tsoil(nlon, nsoilmx, nbsrf)
56 ✗ REAL :: tslab(nlon), seaice(nlon)
57 ✗ REAL fder(nlon)
58
59
60
61 ! Arguments :
62 ! -----------
63
64 ! integer radpas
65 INTEGER it, unit, i, k, itap
66
67 REAL rugos, albedo
68 REAL tsurf
69 REAL time, timestep, day, day0
70 REAL qsol_f
71 ✗ REAL rugsrel(nlon)
72 LOGICAL alb_ocean
73
74 CHARACTER *80 ans, file_forctl, file_fordat, file_start
75 CHARACTER *100 file, var
76 CHARACTER *2 cnbl
77
78 ✗ REAL phy_nat(nlon, year_len)
79 ✗ REAL phy_alb(nlon, year_len)
80 ✗ REAL phy_sst(nlon, year_len)
81 ✗ REAL phy_bil(nlon, year_len)
82 ✗ REAL phy_rug(nlon, year_len)
83 ✗ REAL phy_ice(nlon, year_len)
84 ✗ REAL phy_fter(nlon, year_len)
85 ✗ REAL phy_foce(nlon, year_len)
86 ✗ REAL phy_fsic(nlon, year_len)
87 ✗ REAL phy_flic(nlon, year_len)
88
89 INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
90 !$OMP THREADPRIVATE(read_climoz)
91
92 ! -------------------------------------------------------------------------
93 ! declaration pour l'appel a phyredem
94 ! -------------------------------------------------------------------------
95
96 ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
97 ✗ REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
98 ! real pbl_tke(nlon,llm,nbsrf)
99 ! real rain_fall(nlon),snow_fall(nlon)
100 ! real solsw(nlon), sollw(nlon),radsol(nlon)
101 ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
102 ! real ratqs(nlon,llm)
103 ! real clwcon(nlon,llm)
104
105 INTEGER longcles
106 PARAMETER (longcles=20)
107 REAL clesphy0(longcles)
108
109
110 ! -----------------------------------------------------------------------
111 ! dynamial tendencies :
112 ! ---------------------
113
114 INTEGER l, ierr, aslun
115
116 REAL paire
117
118 ! Local
119 CHARACTER (LEN=20) :: modname='phyaqua'
120 CHARACTER (LEN=80) :: abort_message
121
122
123 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
124 ! INITIALISATIONS
125 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
126
127 ! -----------------------------------------------------------------------
128 ! Initialisations des constantes
129 ! -------------------------------
130
131 !IF (calend .EQ. "earth_360d") Then
132 ✗ year_len_phy = year_len
133 !END IF
134
135 ✗ if (year_len.ne.360) then
136 ✗ write (*,*) year_len
137 ✗ write (*,*) 'iniaqua: 360 day calendar is required !'
138 ✗ stop
139 endif
140
141 ✗ type_aqua = iflag_phys/100
142 ✗ type_profil = iflag_phys - type_aqua*100
143 ✗ PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
144
145 ✗ IF (klon/=nlon) THEN
146 ✗ WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
147 ✗ abort_message= 'probleme de dimensions dans iniaqua'
148 ✗ CALL abort_physic(modname,abort_message,1)
149 END IF
150 ✗ CALL phys_state_var_init(read_climoz)
151
152
153 ✗ read_climoz = 0
154 day0 = 217.
155 day = day0
156 it = 0
157 time = 0.
158
159 ! -----------------------------------------------------------------------
160 ! initialisations de la physique
161 ! -----------------------------------------------------------------------
162
163 ✗ day_ini = day_ref
164 ✗ day_end = day_ini + ndays
165
166 ✗ nbapp_rad = 24
167 ✗ CALL getin_p('nbapp_rad', nbapp_rad)
168
169 ! ---------------------------------------------------------------------
170 ! Creation des conditions aux limites:
171 ! ------------------------------------
172 ! Initialisations des constantes
173 ! Ajouter les manquants dans planete.def... (albedo etc)
174 ✗ co2_ppm = 348.
175 ✗ CALL getin_p('co2_ppm', co2_ppm)
176
177 ✗ solaire = 1365.
178 ✗ CALL getin_p('solaire', solaire)
179
180 ! CALL getin('albedo',albedo) ! albedo is set below, depending on
181 ! type_aqua
182 ✗ alb_ocean = .TRUE.
183 ✗ CALL getin_p('alb_ocean', alb_ocean)
184
185 ✗ WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
186 ✗ WRITE (*, *) 'iniaqua: solaire=', solaire
187 ✗ WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
188
189 ✗ radsol = 0.
190 qsol_f = 10.
191
192 ! Conditions aux limites:
193 ! -----------------------
194
195 ✗ qsol(:) = qsol_f
196 ✗ rugsrel = 0.0 ! (rugsrel = rugoro)
197 ✗ rugoro = 0.0
198 ✗ u_ancien = 0.0
199 ✗ v_ancien = 0.0
200 ✗ agesno = 50.0
201 ! Relief plat
202 ✗ zmea = 0.
203 ✗ zstd = 0.
204 ✗ zsig = 0.
205 ✗ zgam = 0.
206 ✗ zthe = 0.
207 ✗ zpic = 0.
208 ✗ zval = 0.
209
210 ! Une seule surface
211 ✗ pctsrf = 0.
212 ✗ IF (type_aqua==1) THEN
213 ✗ rugos = 1.E-4
214 albedo = 0.19
215 ✗ pctsrf(:, is_oce) = 1.
216 ✗ ELSE IF (type_aqua==2) THEN
217 ✗ rugos = 0.03
218 albedo = 0.1
219 ✗ pctsrf(:, is_ter) = 1.
220 END IF
221
222 ✗ CALL getin_p('rugos', rugos)
223
224 ✗ WRITE (*, *) 'iniaqua: rugos=', rugos
225 ✗ zmasq(:) = pctsrf(:, is_ter)
226
227 ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
228 ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
229
230 ! Si alb_ocean on calcule un albedo oceanique moyen
231 ! if (alb_ocean) then
232 ! Voir pourquoi on avait ca.
233 ! CALL ini_alb_oce(phy_alb)
234 ! else
235 ✗ phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
236 ! endif !alb_ocean
237
238 ✗ DO i = 1, year_len
239 ! IM Terraplanete phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
240 ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
241
242 ✗ phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
243 ✗ phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
244 ✗ phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
245 ✗ phy_ice(:, i) = 0.0 ! fraction de glace (?)
246 ✗ phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
247 ✗ phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
248 ✗ phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
249 ✗ phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
250 END DO
251 ! IM calcul profil sst
252 ✗ CALL profil_sst(nlon, latitude, type_profil, phy_sst)
253
254 ✗ IF (grid_type==unstructured) THEN
255 CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
256 phy_fter, phy_foce, phy_flic, phy_fsic)
257 ELSE
258
259 CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
260 ✗ phy_fter, phy_foce, phy_flic, phy_fsic)
261 ENDIF
262
263 ! ---------------------------------------------------------------------
264 ! Ecriture de l'etat initial:
265 ! ---------------------------
266
267
268 ! Ecriture etat initial physique
269
270 ✗ timestep = pdtphys
271 ✗ radpas = nint(rday/timestep/float(nbapp_rad))
272
273 DO i = 1, longcles
274 ✗ clesphy0(i) = 0.
275 END DO
276 clesphy0(1) = float(iflag_con)
277 clesphy0(2) = float(nbapp_rad)
278 ! IF( cycle_diurne ) clesphy0(3) = 1.
279 clesphy0(3) = 1. ! cycle_diurne
280 clesphy0(4) = 1. ! soil_model
281 clesphy0(5) = 1. ! new_oliq
282 clesphy0(6) = 0. ! ok_orodr
283 clesphy0(7) = 0. ! ok_orolf
284 clesphy0(8) = 0. ! ok_limitvrai
285
286
287 ! =======================================================================
288 ! Profils initiaux
289 ! =======================================================================
290
291 ! On initialise les temperatures de surfaces comme les sst
292 ✗ DO i = 1, nlon
293 ✗ ftsol(i, :) = phy_sst(i, 1)
294 ✗ tsoil(i, :, :) = phy_sst(i, 1)
295 ✗ tslab(i) = phy_sst(i, 1)
296 END DO
297
298 ✗ falbe(:, :) = albedo
299 ✗ falblw(:, :) = albedo
300 ✗ rain_fall(:) = 0.
301 ✗ snow_fall(:) = 0.
302 ✗ solsw(:) = 0.
303 ✗ sollw(:) = 0.
304 ✗ radsol(:) = 0.
305
306 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
307 ! intialisation bidon mais pas grave
308 ✗ t_ancien(:, :) = 0.
309 ✗ q_ancien(:, :) = 0.
310 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
311 ✗ rnebcon = 0.
312 ✗ ratqs = 0.
313 ✗ clwcon = 0.
314 ✗ pbl_tke = 1.E-8
315
316 ! variables supplementaires pour appel a plb_surface_init
317 ✗ fder(:) = 0.
318 ✗ seaice(:) = 0.
319 ✗ run_off_lic_0 = 0.
320 ✗ fevap = 0.
321
322
323 ! Initialisations necessaires avant phyredem
324 ✗ type_ocean = 'force'
325 ✗ CALL fonte_neige_init(run_off_lic_0)
326 ✗ qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
327 ✗ snsrf(:, :) = 0. ! couverture de neige des sous surface
328 ✗ z0m(:, :) = rugos ! couverture de neige des sous surface
329 ✗ z0h=z0m
330
331
332 ✗ CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
333
334 ✗ PRINT *, 'iniaqua: before phyredem'
335
336 ✗ pbl_tke(:,:,:) = 1.e-8
337 ✗ falb1 = albedo
338 ✗ falb2 = albedo
339 ✗ zmax0 = 0.
340 ✗ f0 = 0.
341 ✗ sig1 = 0.
342 ✗ w01 = 0.
343 ✗ wake_deltat = 0.
344 ✗ wake_deltaq = 0.
345 ✗ wake_s = 0.
346 ✗ wake_dens = 0.
347 ✗ wake_cstar = 0.
348 ✗ wake_pe = 0.
349 ✗ wake_fip = 0.
350 ✗ fm_therm = 0.
351 ✗ entr_therm = 0.
352 ✗ detr_therm = 0.
353 ✗ ale_bl = 0.
354 ✗ ale_bl_trig =0.
355 ✗ alp_bl =0.
356 ✗ treedrg(:,:,:)=0.
357
358 ✗ u10m = 0.
359 ✗ v10m = 0.
360
361 ✗ ql_ancien = 0.
362 ✗ qs_ancien = 0.
363 ✗ u_ancien = 0.
364 ✗ v_ancien = 0.
365 ✗ prw_ancien = 0.
366 ✗ prlw_ancien = 0.
367 ✗ prsw_ancien = 0.
368
369 ✗ ale_wake = 0.
370 ✗ ale_bl_stat = 0.
371
372
373 !ym error : the sub surface dimension is the third not second : forgotten for iniaqua
374 ! falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
375 ! falb_dir(:,is_oce,:)=0.5; falb_dir(:,is_sic,:)=0.6
376 ✗ falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
377 ✗ falb_dir(:,:,is_oce)=0.5; falb_dir(:,:,is_sic)=0.6
378
379 !ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
380 !ym probably the uninitialized value was 0 for standard (regular grid) case
381 ✗ falb_dif(:,:,:)=0
382
383
384 ✗ CALL phyredem('startphy.nc')
385
386 ✗ PRINT *, 'iniaqua: after phyredem'
387 ✗ CALL phys_state_var_end
388
389 ✗ RETURN
390 END SUBROUTINE iniaqua
391
392
393 ! ====================================================================
394 ! ====================================================================
395 ✗ SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
396 ✗ USE dimphy
397 IMPLICIT NONE
398 ! ====================================================================
399 ! =============================================================
400 ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
401 ! Auteur : A. Campoy et F. Hourdin
402 ! Objet : calculer les valeurs moyennes du cos de l'angle zenithal
403 ! et l'ensoleillement moyen entre gmtime1 et gmtime2
404 ! connaissant la declinaison, la latitude et la longitude.
405
406 ! Dans cette version particuliere, on calcule le rayonnement
407 ! moyen sur l'année à chaque latitude.
408 ! angle zenithal calculé pour obtenir un
409 ! Fit polynomial de l'ensoleillement moyen au sommet de l'atmosphere
410 ! en moyenne annuelle.
411 ! Spécifique de la terre. Utilisé pour les aqua planetes.
412
413 ! Rque : Different de la routine angle en ce sens que zenang
414 ! fournit des moyennes de pmu0 et non des valeurs
415 ! instantanees, du coup frac prend toutes les valeurs
416 ! entre 0 et 1.
417 ! Date : premiere version le 13 decembre 1994
418 ! revu pour GCM le 30 septembre 1996
419 ! ===============================================================
420 ! longi----INPUT : la longitude vraie de la terre dans son plan
421 ! solaire a partir de l'equinoxe de printemps (degre)
422 ! gmtime---INPUT : temps universel en fraction de jour
423 ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
424 ! lat------INPUT : latitude en degres
425 ! long-----INPUT : longitude en degres
426 ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
427 ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
428 ! ================================================================
429 include "YOMCST.h"
430 ! ================================================================
431 LOGICAL cycle_diurne
432 REAL gmtime
433 REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
434 ! ================================================================
435 INTEGER i
436 REAL gmtime1, gmtime2
437 REAL pi_local
438
439
440 ✗ REAL rmu0m(klon), rmu0a(klon)
441
442
443 pi_local = 4.0*atan(1.0)
444
445 ! ================================================================
446 ! Calcul de l'angle zenithal moyen sur la journee
447 ! ================================================================
448
449 ✗ DO i = 1, klon
450 ✗ fract(i) = 1.
451 ! Calcule du flux moyen
452 ✗ IF (abs(rlat(i))<=28.75) THEN
453 ✗ rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
454 ✗ ELSE IF (abs(rlat(i))<=43.75) THEN
455 ✗ rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
456 ✗ ELSE IF (abs(rlat(i))<=71.25) THEN
457 ✗ rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
458 ELSE
459 ✗ rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
460 END IF
461 END DO
462
463 ! ================================================================
464 ! Avec ou sans cycle diurne
465 ! ================================================================
466
467 ✗ IF (cycle_diurne) THEN
468
469 ! On redecompose flux au sommet suivant un cycle diurne idealise
470 ! identique a toutes les latitudes.
471
472 ✗ DO i = 1, klon
473 ✗ rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
474 rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
475 ✗ rmu0a(i)/sqrt(2.)
476 END DO
477
478 ✗ DO i = 1, klon
479 ✗ IF (rmu0(i)<=0.) THEN
480 ✗ rmu0(i) = 0.
481 ✗ fract(i) = 0.
482 ELSE
483 ✗ fract(i) = 1.
484 END IF
485 END DO
486
487 ! Affichage de l'angel zenitale
488 ! print*,'************************************'
489 ! print*,'************************************'
490 ! print*,'************************************'
491 ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
492 ! print*,'rmu0m=',rmu0m(i)
493 ! print*,'rmu0a=',rmu0a(i)
494 ! print*,'rmu0=',rmu0(i)
495
496 ELSE
497
498 ✗ DO i = 1, klon
499 ✗ fract(i) = 0.5
500 ✗ rmu0(i) = rmu0m(i)*2.
501 END DO
502
503 END IF
504
505 ✗ RETURN
506 END SUBROUTINE zenang_an
507
508 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
509
510 ✗ SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
511 phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
512
513 USE mod_phys_lmdz_para, ONLY: is_omp_master, klon_mpi
514 USE mod_phys_lmdz_transfert_para, ONLY: gather_omp
515 IMPLICIT NONE
516
517 include "netcdf.inc"
518
519 INTEGER, INTENT (IN) :: klon
520 REAL, INTENT (IN) :: phy_nat(klon, 360)
521 REAL, INTENT (IN) :: phy_alb(klon, 360)
522 REAL, INTENT (IN) :: phy_sst(klon, 360)
523 REAL, INTENT (IN) :: phy_bil(klon, 360)
524 REAL, INTENT (IN) :: phy_rug(klon, 360)
525 REAL, INTENT (IN) :: phy_ice(klon, 360)
526 REAL, INTENT (IN) :: phy_fter(klon, 360)
527 REAL, INTENT (IN) :: phy_foce(klon, 360)
528 REAL, INTENT (IN) :: phy_flic(klon, 360)
529 REAL, INTENT (IN) :: phy_fsic(klon, 360)
530
531 REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
532 ! on the whole physics grid
533
534 ✗ END SUBROUTINE writelim_unstruct
535
536
537
538 ✗ SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
539 phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
540
541 USE mod_phys_lmdz_para, ONLY: is_master
542 USE mod_grid_phy_lmdz, ONLY: klon_glo
543 USE mod_phys_lmdz_transfert_para, ONLY: gather
544 USE phys_cal_mod, ONLY: year_len
545 IMPLICIT NONE
546 include "netcdf.inc"
547
548 INTEGER, INTENT (IN) :: klon
549 REAL, INTENT (IN) :: phy_nat(klon, year_len)
550 REAL, INTENT (IN) :: phy_alb(klon, year_len)
551 REAL, INTENT (IN) :: phy_sst(klon, year_len)
552 REAL, INTENT (IN) :: phy_bil(klon, year_len)
553 REAL, INTENT (IN) :: phy_rug(klon, year_len)
554 REAL, INTENT (IN) :: phy_ice(klon, year_len)
555 REAL, INTENT (IN) :: phy_fter(klon, year_len)
556 REAL, INTENT (IN) :: phy_foce(klon, year_len)
557 REAL, INTENT (IN) :: phy_flic(klon, year_len)
558 REAL, INTENT (IN) :: phy_fsic(klon, year_len)
559
560 ✗ REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
561 ! on the whole physics grid
562 INTEGER :: k
563 INTEGER ierr
564 INTEGER dimfirst(3)
565 INTEGER dimlast(3)
566
567 INTEGER nid, ndim, ntim
568 INTEGER dims(2), debut(2), epais(2)
569 INTEGER id_tim
570 INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
571 INTEGER id_fter, id_foce, id_fsic, id_flic
572
573 ✗ IF (is_master) THEN
574
575 ✗ PRINT *, 'writelim: Ecriture du fichier limit'
576
577 ✗ ierr = nf_create('limit.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
578
579 ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
580 ✗ 'Fichier conditions aux limites')
581 ! ! ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
582 ✗ ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
583 ✗ ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
584
585 ✗ dims(1) = ndim
586 ✗ dims(2) = ntim
587
588 ✗ ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
589 ✗ ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
590
591 ✗ ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
592 ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
593 ✗ 'Nature du sol (0,1,2,3)')
594
595 ✗ ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
596 ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
597 ✗ 'Temperature superficielle de la mer')
598
599 ✗ ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
600 ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
601 ✗ 'Reference flux de chaleur au sol')
602
603 ✗ ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
604 ✗ ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
605
606 ✗ ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
607 ✗ ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
608
609 ✗ ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
610 ✗ ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
611 ✗ ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
612 ✗ ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
613 ✗ ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
614 ✗ ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
615 ✗ ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
616 ✗ ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
617
618 ✗ ierr = nf_enddef(nid)
619 ✗ IF (ierr/=nf_noerr) THEN
620 ✗ WRITE (*, *) 'writelim error: failed to end define mode'
621 ✗ WRITE (*, *) nf_strerror(ierr)
622 END IF
623
624
625 ! write the 'times'
626 ✗ DO k = 1, year_len
627 ✗ ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
628 ✗ IF (ierr/=nf_noerr) THEN
629 ✗ WRITE (*, *) 'writelim error with temps(k),k=', k
630 ✗ WRITE (*, *) nf_strerror(ierr)
631 END IF
632 END DO
633
634 END IF ! of if (is_master)
635
636 ! write the fields, after having collected them on master
637
638 ✗ CALL gather(phy_nat, phy_glo)
639 ✗ IF (is_master) THEN
640 ✗ ierr = nf_put_var_double(nid, id_nat, phy_glo)
641 ✗ IF (ierr/=nf_noerr) THEN
642 ✗ WRITE (*, *) 'writelim error with phy_nat'
643 ✗ WRITE (*, *) nf_strerror(ierr)
644 END IF
645 END IF
646
647 ✗ CALL gather(phy_sst, phy_glo)
648 ✗ IF (is_master) THEN
649 ✗ ierr = nf_put_var_double(nid, id_sst, phy_glo)
650 ✗ IF (ierr/=nf_noerr) THEN
651 ✗ WRITE (*, *) 'writelim error with phy_sst'
652 ✗ WRITE (*, *) nf_strerror(ierr)
653 END IF
654 END IF
655
656 ✗ CALL gather(phy_bil, phy_glo)
657 ✗ IF (is_master) THEN
658 ✗ ierr = nf_put_var_double(nid, id_bils, phy_glo)
659 ✗ IF (ierr/=nf_noerr) THEN
660 ✗ WRITE (*, *) 'writelim error with phy_bil'
661 ✗ WRITE (*, *) nf_strerror(ierr)
662 END IF
663 END IF
664
665 ✗ CALL gather(phy_alb, phy_glo)
666 ✗ IF (is_master) THEN
667 ✗ ierr = nf_put_var_double(nid, id_alb, phy_glo)
668 ✗ IF (ierr/=nf_noerr) THEN
669 ✗ WRITE (*, *) 'writelim error with phy_alb'
670 ✗ WRITE (*, *) nf_strerror(ierr)
671 END IF
672 END IF
673
674 ✗ CALL gather(phy_rug, phy_glo)
675 ✗ IF (is_master) THEN
676 ✗ ierr = nf_put_var_double(nid, id_rug, phy_glo)
677 ✗ IF (ierr/=nf_noerr) THEN
678 ✗ WRITE (*, *) 'writelim error with phy_rug'
679 ✗ WRITE (*, *) nf_strerror(ierr)
680 END IF
681 END IF
682
683 ✗ CALL gather(phy_fter, phy_glo)
684 ✗ IF (is_master) THEN
685 ✗ ierr = nf_put_var_double(nid, id_fter, phy_glo)
686 ✗ IF (ierr/=nf_noerr) THEN
687 ✗ WRITE (*, *) 'writelim error with phy_fter'
688 ✗ WRITE (*, *) nf_strerror(ierr)
689 END IF
690 END IF
691
692 ✗ CALL gather(phy_foce, phy_glo)
693 ✗ IF (is_master) THEN
694 ✗ ierr = nf_put_var_double(nid, id_foce, phy_glo)
695 ✗ IF (ierr/=nf_noerr) THEN
696 ✗ WRITE (*, *) 'writelim error with phy_foce'
697 ✗ WRITE (*, *) nf_strerror(ierr)
698 END IF
699 END IF
700
701 ✗ CALL gather(phy_fsic, phy_glo)
702 ✗ IF (is_master) THEN
703 ✗ ierr = nf_put_var_double(nid, id_fsic, phy_glo)
704 ✗ IF (ierr/=nf_noerr) THEN
705 ✗ WRITE (*, *) 'writelim error with phy_fsic'
706 ✗ WRITE (*, *) nf_strerror(ierr)
707 END IF
708 END IF
709
710 ✗ CALL gather(phy_flic, phy_glo)
711 ✗ IF (is_master) THEN
712 ✗ ierr = nf_put_var_double(nid, id_flic, phy_glo)
713 ✗ IF (ierr/=nf_noerr) THEN
714 ✗ WRITE (*, *) 'writelim error with phy_flic'
715 ✗ WRITE (*, *) nf_strerror(ierr)
716 END IF
717 END IF
718
719 ! close file:
720 ✗ IF (is_master) THEN
721 ✗ ierr = nf_close(nid)
722 END IF
723
724 ✗ END SUBROUTINE writelim
725
726 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
727
728 ✗ SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
729 USE dimphy
730 USE phys_cal_mod , ONLY: year_len
731 IMPLICIT NONE
732
733 INTEGER nlon, type_profil, i, k, j
734 REAL :: rlatd(nlon), phy_sst(nlon, year_len)
735 INTEGER imn, imx, amn, amx, kmn, kmx
736 INTEGER p, pplus, nlat_max
737 PARAMETER (nlat_max=72)
738 REAL x_anom_sst(nlat_max)
739 CHARACTER (LEN=20) :: modname='profil_sst'
740 CHARACTER (LEN=80) :: abort_message
741
742 ✗ IF (klon/=nlon) THEN
743 ✗ abort_message='probleme de dimensions dans profil_sst'
744 ✗ CALL abort_physic(modname,abort_message,1)
745 ENDIF
746 ✗ WRITE (*, *) ' profil_sst: type_profil=', type_profil
747 ✗ DO i = 1, year_len
748 ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
749
750 ! Rajout fbrlmd
751
752 ✗ IF (type_profil==1) THEN
753 ! Méthode 1 "Control" faible plateau à l'Equateur
754 ✗ DO j = 1, klon
755 ✗ phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
756 ! PI/3=1.047197551
757 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
758 ✗ phy_sst(j, i) = 273.
759 END IF
760 END DO
761 END IF
762 ✗ IF (type_profil==2) THEN
763 ! Méthode 2 "Flat" fort plateau à l'Equateur
764 ✗ DO j = 1, klon
765 ✗ phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
766 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
767 ✗ phy_sst(j, i) = 273.
768 END IF
769 END DO
770 END IF
771
772
773 ✗ IF (type_profil==3) THEN
774 ! Méthode 3 "Qobs" plateau réel à l'Equateur
775 ✗ DO j = 1, klon
776 phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
777 ✗ rlatd(j))**4)
778 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
779 ✗ phy_sst(j, i) = 273.
780 END IF
781 END DO
782 END IF
783
784 ✗ IF (type_profil==4) THEN
785 ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
786 ✗ DO j = 1, klon
787 phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
788 ✗ rlatd(j))**4)
789 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
790 ✗ phy_sst(j, i) = 273.
791 END IF
792 END DO
793 END IF
794
795 ✗ IF (type_profil==5) THEN
796 ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
797 ✗ DO j = 1, klon
798 phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
799 ✗ *rlatd(j))**4)
800 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
801 ✗ phy_sst(j, i) = 273. + 2.
802 END IF
803
804 END DO
805 END IF
806
807 ✗ IF (type_profil==6) THEN
808 ! Méthode 6 "cst" valeur constante de SST
809 ✗ DO j = 1, klon
810 ✗ phy_sst(j, i) = 288.
811 END DO
812 END IF
813
814
815 ✗ IF (type_profil==7) THEN
816 ! Méthode 7 "cst" valeur constante de SST +2
817 ✗ DO j = 1, klon
818 ✗ phy_sst(j, i) = 288. + 2.
819 END DO
820 END IF
821
822 p = 0
823 ✗ IF (type_profil==8) THEN
824 ! Méthode 8 profil anomalies SST du modèle couplé AR4
825 ✗ DO j = 1, klon
826 ✗ IF (rlatd(j)==rlatd(j-1)) THEN
827 phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
828 ✗ 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
829 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
830 ✗ phy_sst(j, i) = 273. + x_anom_sst(pplus)
831 END IF
832 ELSE
833 ✗ p = p + 1
834 ✗ pplus = 73 - p
835 phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
836 ✗ 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
837 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
838 ✗ phy_sst(j, i) = 273. + x_anom_sst(pplus)
839 END IF
840 ✗ WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
841 END IF
842 END DO
843 END IF
844
845 ✗ IF (type_profil==9) THEN
846 ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
847 ✗ DO j = 1, klon
848 phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
849 ✗ *rlatd(j))**4)
850 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
851 ✗ phy_sst(j, i) = 273. - 2.
852 END IF
853 END DO
854 END IF
855
856
857 ✗ IF (type_profil==10) THEN
858 ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
859 ✗ DO j = 1, klon
860 phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
861 ✗ *rlatd(j))**4)
862 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
863 ✗ phy_sst(j, i) = 273. + 4.
864 END IF
865 END DO
866 END IF
867
868 ✗ IF (type_profil==11) THEN
869 ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
870 ✗ DO j = 1, klon
871 phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
872 ✗ rlatd(j))**4)
873 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
874 ✗ phy_sst(j, i) = 273.
875 END IF
876 END DO
877 END IF
878
879 ✗ IF (type_profil==12) THEN
880 ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
881 ✗ DO j = 1, klon
882 phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
883 ✗ *rlatd(j))**4)
884 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
885 ✗ phy_sst(j, i) = 273. + 4.
886 END IF
887 END DO
888 END IF
889
890 ✗ IF (type_profil==13) THEN
891 ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
892 ✗ DO j = 1, klon
893 phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
894 ✗ rlatd(j))**4)
895 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
896 ✗ phy_sst(j, i) = 273.
897 END IF
898 END DO
899 END IF
900
901 ✗ IF (type_profil==14) THEN
902 ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
903 ✗ DO j = 1, klon
904 phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
905 ✗ *rlatd(j))**4)
906 ✗ IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
907 ✗ phy_sst(j, i) = 273.
908 END IF
909 END DO
910 END IF
911
912 ✗ if (type_profil.EQ.20) then
913 ✗ print*,'Profile SST 20'
914 ! Méthode 13 "Qmax2K" plateau réel �| l'Equateur augmenté +2K
915
916 ✗ do j=1,klon
917 ✗ phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
918 enddo
919 endif
920
921 ✗ if (type_profil.EQ.21) then
922 ✗ print*,'Profile SST 21'
923 ! Méthode 13 "Qmax2K" plateau réel �| l'Equateur augmenté +2K
924 ✗ do j=1,klon
925 ✗ phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
926 enddo
927 endif
928
929
930
931 END DO
932
933 ! IM beg : verif profil SST: phy_sst
934 ✗ amn = min(phy_sst(1,1), 1000.)
935 ✗ amx = max(phy_sst(1,1), -1000.)
936 ✗ imn = 1
937 ✗ kmn = 1
938 ✗ imx = 1
939 ✗ kmx = 1
940 ✗ DO k = 1, year_len
941 ✗ DO i = 2, nlon
942 ✗ IF (phy_sst(i,k)<amn) THEN
943 ✗ amn = phy_sst(i, k)
944 ✗ imn = i
945 ✗ kmn = k
946 END IF
947 ✗ IF (phy_sst(i,k)>amx) THEN
948 ✗ amx = phy_sst(i, k)
949 ✗ imx = i
950 ✗ kmx = k
951 END IF
952 END DO
953 END DO
954
955 ✗ PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
956 ✗ PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
957 ! IM end : verif profil SST: phy_sst
958
959 ✗ RETURN
960 END SUBROUTINE profil_sst
961
962 END MODULE phyaqua_mod
963