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