source: LMDZ5/trunk/libf/phylmd/phyaqua_mod.F90 @ 5451

Last change on this file since 5451 was 2979, checked in by Ehouarn Millour, 7 years ago

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